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PARAMLTRIC  SCALING  LAWS 

) 

Part  II:  A Computer  Model  for  Deriving  the  Performance 

Characteristics  of  Saturation-Limited  Parametric  Sonars 

i by 

F.  H.  Fenlon  and  J.  W.  Kesner 


SUMMARY 

The  purpose  and  significance  of  this  ARPA  sponsored  report  is 
to  provide  a detailed  description  and  listing  of  the  computer  programs 
developed  to  compute  the  amplitude,  frequency,  and  angular  response 
characteristics  of  parametric  sonars,  from  an  analytical  model  developed 
in  Part  I.  Technical  results  include  examples  of  the  performance 
characteristics  and  beam  patterns  generated  by  the  program.  Implications 
for  further  research  would  be  to  work  backwards  from  the  asymptotic  far— 
field  solutions  derived  in  this  report  to  compute  the  range  dependence 
of  the  entire  spectrum  throughout  the  nonlinear  interaction  zone  via  the 


paraxial  wave  equation. 


TABLE  OF  CONTENTS 


LIST  OF  SYMBOLS  

THEORY  

PROGRAM  INPUT  

PROGRAM  OUTPUT  

REFERENCES  

APPENDIX  A 

FIGURE  1 (Pages  37-85)  . . . . 

FIGURE  2 (Pages  87-99)  . . . . 

FIGURE  3 (Pages  101-104)  . . . 

PROGRAM  LISTING  (Pages  105-118) 


. 1< 

. 11 
. 13 
. 37 
. 8] 

. 101 
. 105 


lb 


> 

A 


1. 1ST  OF  SYMBOLS 

density  of  the  medium 
small  signal  speed  of  sound 
i on linen r parameter  for  liquids 
in  gases 

first  and  second-order  coefficients  in  the  equation 
of  state 

specific  heats  at  constant  pressure  and  constant  volume 

Rayleigh  distance,  and  end-fire  truncation  distance 

primary  wave  frequencies 

mean  primary  wave  frequency 

sum  cr  difference  frequency 

nrimary  wave  numbers 

mean  primary  wave  number 

ms  primary  wave  amplitudes  at  the  source 

nondispersive  thermo  viscous  attenuation  parameter 

attenuation  coefficient 

primary  wave  attenuation  coefficients 

sum  or  difference-frequency  attenuation  coefficient 

combined  attenuation  coefficient 

mean  primary  wave  attenuation  coefficient 
mean  primary  wave  acoustic  Reynolds  number 
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E,  (a  r ) e 
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Vo 


El<Vi>  C 


Vo 


1 + 


log 1()(*/f ') 


lo8l0(fo/f-T 


E^z) 


-x 


dx 


D1(0),(i-1,2) 

D+(0) 


exponential  Integral 

far-field  primary  wave  directivity  functions 


far-field  sum  or  difference-frequency  directivity 
function  for  a "spreading-loss-limited"  parametric 
array 


M°) 


far-field  sum  or  difference-frequency  directivity 
function  for  an  "absorption-limited"  parametric  array 


b,(0) 

dt(0) 

DI 

SL 


convolved  beam  pattern 
I),  (0),  or  b+(0) 
directivity  index 

combined  rms  primary  wave  source  level  in  dB//lyPam 


SL 


equivalent  mean  rms  primary  wave  source  level  at  1 m 
in  dB//lyPam 


SL. 


equivalent  sum  or  d if f erence-f r eouency  rms  source 
level  at  1 n in  dB//lyPan 


SL*  = SL  -*-20  log  nf  , (f  in  kHz)  scaled  combined  rms  primary  wave  source 
° 0 10  o o level  in  dB//luPam  kHz 


SL*  « SL  +20  log  _f  scaled  equivalent  mean  rms  primary  wave  source 

00  oo  LG  O - . . in  / il  n l.tt. 


level  at  1 m in  dB//lyPam  kHz 


SL*  * SL  +20  log^QfQ  scaled  equivalent  sum  or  difference  frequency  rms 


source  level  in  dB//lyPa  k!!z 
20  log^Qq+  = SL+  - SLq  parametric  (sum  or  difference)  conversion  efficiency 
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I ll  KORY 

In  Part  I of  this  report  an  analytical  solution  for  the  far-field 
pressure  of  a parametric  array  was  derived  from  Burgers'  equation.  Section  2 
of  Part  i gave  this  solution  in  terms  of  scaled  source  levels.  Since  these 
scaled  expressions  form  the  basis  of  the  computer  program  outlined  in  this 
part  of  the  report,  they  are  reexpressed  here  for  the  benefit  of  the  reader 
as  follows: 

Let  SL  ■ the  combined  orimary  wave  source  level  at  lm,  in  dB//lpPam 
o 

and  SL+  = the  equivalent  sum  or  di f f erence-f reouencv  source  level 


referred  to  lm,  in  dB// IpPam. 

Then  the  scaled  source  levels  SL*  and  SL*  are  defined  as, 

o * 


SL*  • SL  +20  loglnf  ) 

0 0 10  ° 1 f in  kHz 

(la) 

sl;  - sl;  + 20  ioe10fo  J 

(lb) 

where 

f0  ’ 7(tl  + £2>- 

(2) 

f^  and  f ^ being  the  operating  frequencies  of  the  parametric  sonar.  fQ  is 
thus  the  mean  primary  v/ave  frequencv. 

Parametric  Transmitter: 

From  Section  2 of  Part  I the  conversion  efficiencv  q_  of  a 
parametric  transmitter  can  be  expressed  in  logarithmic  form  as, 

20  log10n_  = (SL_  - SLq) 

= (SL*  - SL*),  from  Eqs.  (la)  and  (lb) 

= 20  log1QF  - 2 On  log1()(fo/f_)  (3) 
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,1<X1>I1^X^) 


K - (2/xo)  lo(X])lo(x2)  : *<’ 


20  1o8,0ao  " SLS  + 20  loRl(/  + 20  loR10N 
20  log^Xj  * SL*  + 20  log1Q^+  20  log1()N,  (i-1,2) 


N - (2/2*  B/p  c3)(103) 
o o 


t ■ E1<tITCo'eXp<“Tco)!  r'  ’ F'l*“Tto'<!Xp<°Tro'! 


v;  ■ v.'V1-1 


n = 1 + 


log10(^/^') 


log10(fo/f-V 


The  program  computes,  plots,  and  tabulates  the  conversion  efficiency 


20  logn  as  a function  of  the  scaled  combined  primary  wave  source  level 


SL*  for  fixed  values  of  a_r  and  f /f  . It  also  plots  and  tabulates  the 


scaled  difference-frequency  source  level  SL*  is  a function  of  SL*  for 

fixed  values  of  oi_r  and  f / f . In  addition  the  frequenev  response  index 
To  o - 


n,  as  defined  by  F.q.  8 is  tabulated.  Rv  reading  in  the  relative  amplitude 


n /p  of  the  rms  primary  wave  amplitudes  at  the  source  the  program  will 
ol  o2 


provide  results  for  cases  where  the  drive  amplitudes  are  unequal. 


Parametric  Receiver: 


For  a limited  class  of  problems  such  as  that  of  a plane  wave 


low  frequency  signal  interacting  with  plane  wave-  'absorption  limited 
finite-amplitude  source,  or  a spherical  wave  signal  interacting  with  a 


spherical  nump  wave,  the  same  expressions  used  to  define  the  performance 


of  a parametric  transmitter  will  also  define  the  performance  of  a 


parametric  receiver  if  f Q / f _ is  replaced  by  f Q/f+  in  F.q.  3 and  if 


2 


iv  7 Tin*  roison  for  tills  csn  be  seen  by 
a r is  replaced  by  -aTro  in  bo.  7.  The  reason 


writing  oy^  as  aT+r0  wliere' 


Vo 


(<*1  + a2  “ Vc 


6 { f ^ + f2  _ ^fl  ~ ^r<: 


- ? 26flf2ro- 

Thus  cT  ro  - 26f1f2ro,  for  the  difference-frequency  signal  (10a) 


*T  ro  "'2^flf2ro 
+ 


-a  r . for  the  sum-frequency  signal, 
T o 


(10b) 


In  the  computer  output  20  log1Qn_  becomes  20  log1Qn+  and  SL_  becomes  + 

The  substlcucton  of  x - -Vo  in  the  exmeentlal  integral  EjU)  transforms 

1 

this  to  Ei(z),  as  defined  bv  Abramowitz  and  Segun. 

It  should  be  noted  :Ln  the  ease  of  a narametrlc  receiver  that 
the  index  n is  still  calculate  by  in-  8 since  the  sun,  frequency  component 
must  have  the  same  frequency  response  as  the  difference-frequency  component. 


Evaluation  of  E1(z)cxp(z)  fov  z >,  10.0: 

For  z > 10.0,  the  function  E1(z)exp(z)  which  is  used  in 
calculating  l,  i\  and  lQ  as  .iefined  by  Eqs.  7 and  22,  respectively,  ls^ 
determined  from  an  approximate  expression  given  b-  Ahramowitr  and  Segun  as, 

z2  + A^z  + 

E.  (z)exp(z)  = — ^ 5 

^ z + B^z  + B 2Z 


T 


where 


* 4.03640 


- 1.15198 
Bj  - 5.03637 
B2  - 4.19160 


Far-Fleld  Directivity  Function: 


,-2 


For  a "spread ing-loss- 1 imi ted"  parametric  array  (u^r^  .<  10  dB) 
the  program  gives  the  far-field  directivity  function  D4(0)  as, 

,lllx1°1(">l'1lx2D2(0)l,/,I1(x1)l1fa2), 

!<  ’ 1o'Xid1<»)1,„Ix2D2<0)'/  WW  ’ 

where  D^(o),  (1*1,2)  are  the  far-field  primary  wave  directivity  functions 
given  in  the  program  as, 

2.1  (k  a slnO) 

D (0)  * — r — — 7— ; (i*l,2),  for  a circular  piston  (12) 

x k.3  sinn  r t > 

i of  radius  a 


sin(kja  sinO) 
k a sinn 


; (i— 1,2),  for  a line  arrav  of 
length  2a. 


(13) 


_2 

It  should  be  noted  that  although  Eo.  11  is  onlv  applicable  for  a^r^  ^ 10 

the  program  permits  it  to  be  evaluated  for  all  values  of  ot_r'.  For  values 

i o 

_2 

of  ajr^  > 10  dR  therefore,  the  "absorption-limited,"  "virtual-cnd-f ire- 
array"  directivity  function  5,(0)  should  be  used  instead.  This  function 
is  given  in  the  program  bv  the  expression. 


5,(0)  - { y 


a 


n 


n=o  (n+l)a^  rQ  + j2k,ro  sin^(0/2)/  n=o  ^n+^lT  ro 

4 

sin((n+l)tan  *(4/T  )] 

where  a » (T_/r)n  — 

n o / n+1 

|/l  + (ro/4)^i 


} 


(14) 


(15) 


4 


and  20  log10ro  - SI.*  - 20  l.T|0l  y~‘'l  * 


ol 


02  ’ 


(16) 


It  should  be  noted  that  unlike  Eq.  11,  Eq.  14  only  holds  for 

2 

equal  primary  wave  amplitudes.  Although  it  resembles  Bartram's  end-fire- 

directivitv  function,  Eq.  14  is  a new  and  more  general  form  of  the 

3 

directivity  function  described  bv  Childs.  It  was  derived  from  Kuznetsov's 

equation^  via  Burgers'  equation  and  will  bo  the  subject  of  a later  article. 
_2 

For  10  *:  aTro  ' t*ie  conv°lution  °1  Eqs.  H ar>d  14,  which 

is  available  in  the  program,  should  be  used  to  obtain  the  far-field  para- 
metric array  beam  pattern  bf  (9 ) where, 

fir/ 2 


b+  (9) 


D (9  ' )5  + (9 - 9 ' )cos9  ' d9  ' 


(17) 


J-n/2  " 

where  D^(0)  is  defined  by  Eq.  11  and  <$  + (9)  is  given  in  Eq.  14. 
Directivity  Index: 

The  directivity  index  DI  is  calculated  in  the  program  as  a 

function  of  the  scaled  source  level  SL*  for  fixed  values  of  a_r  and 

o To 

fQ/f_,  from  the  expression. 


DI 


10  lo810 


2 

f | d+  (9  ' )d*  (o  ' ) | sinO  ’ dfl' 


(18) 


1 o 

where  d+(9)  stands  for  D+(0),  <5^(0),  or  bf(9)  whichever  is  appropriate 
and  d*(9)  is  its  complex  conjugate. 


Far-Fleld  Monofrequency  Source  Level : 

The  program  can  also  be  used  to  obtain  the  referred  far-field 

scaled  source  level  at  1 m,  SL*  of  a monof rcquencv  source  of  frequency 

f , as  a function  of  the  scaled  "inDut"  source  level  SI.*,  for  fixed 
o’  o’ 

values  of  a r , 
o o 


5 


where 


(19) 


Sl>  « SL*  + 20  1 l()l  * 

'l(x0) 

F'  ■ <2/*°>W 

with  20  log10XQ  - SL*  + 20  loSl0/o  + 20  log1()N 

*o  “ El(2Vo)‘‘Xp(2Vo)- 

F,qs.  19  to  22  will  be  the  subject  of  a later  article.  Note  tha 
ol*  - SL  + 20  log.„f  ; (f  in  kHz), 

i or  10  O O 

where  SL  will  differ  from  SL  as  the  latter  increases,  due  to 

m O 

of  energy  from  f into  self-generated  harmonics  in  the  medium, 
case  of  a parametric  source  the  program  gives  SL*  vs  SL*  as  an 
output"  relationship  for  the  mean  primary  wave  frequency. 


(20) 

(21) 

(22) 

♦* 

(23) 

the  transfer 
In  the 
"input- 
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PROCRAM  INPUT 

The  data  for  the  program  is  in  free  field  format  which  moans 
that  the  data  is  not  restricted  to  anv  particular  columns  of  the  data 
card.  Input  variables  which  begin  with  the  letters  I,  J,  K,  L,  M,  N 
are  integers  (no  decimal  point).  All  others  are  real  variables  and  must 
have  a decimal  point.  The  data  values  on  a card  must  be  separated  from 
each  other  bv  commas. 

Card  1 :I0C  * the  number  of  separate  cases  being  run. 

Card  2 IOP,  IOE,  IOS,  IOD,  100 

Kncli  of  these  integer  indicators  must  be  either  0 or  1 . 

A zero  indicates  that  that  particular  part  of  the  program 
can  be  skipped.  A one  indicates  that  that  section  is  to 
be  done.  IOP  is  for  the  pattern  section.  IOE  is  for  the 
20  log  (ETA)  section.  IOS  is  for  the  SPL-*  section. 

10D  is  for  the  directivitv  index  section.  100  is  for  the 
SPL  (OUT)  section. 

Card  3 ALTRO,  P1P2 

ALTRO  * 0ljro  dB.  ALTRO  is  positive  for  difference- 
freauenev  cases  and  is  negative  for  sum  frecuenev  cases. 
P1P2  is  the  ratio  (P1/P2)  of  the  primary  wave  pressure 
implitudes  of  the  two  drive  frequencies  FI  and  F2  at  the 
source. 
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Card  4 NFM,  NFMRD 


NFM  * the  nunber  of  FO/F-  values.  NFM  must  be  5 or  less. 

NFMKI)  * 0 or  1.  A zero  will  generate  a set  of  FO/F- 

K-l 

cases  accord inp,  to  the  formula  FO/F-  * 5.0*2  , K*l,...,  NFM. 

If  NFMRD  » 1 then  NFM  values  for  FO/F-  must  be  read  in. 

Card  5 FFM(K) , K - 1,  . . . , NFM 

FFM(K)  - the  K1'1  value  of  F0/F-.  if  NFMRD  - 0 this  card 
must  be  omitted. 

Card  6 LIN,  IOC 

LIN  - 0 for  a piston,  ■ l for  a rectangle. 

IOC  “ 1 for  the  "1"  Bessel  function  formula  for  the 
beam  pattern. 

* 2 for  the  beam  pattern  of  the  delta  formula. 

- 3 for  the  beam  pattern  of  the  convolution  of  the 
delta  formula  with  the  "i"  Bessel  function  formula.  If 
1 0 1»  » 0 and  1 01)  = 0 then  tills  card  must  be  omitted. 

Card  7 AKOR  if  LIN  - 0 

AKOL,  AKOW  if  LIN  - 1 

AKOR  - kQa  - 2.0  * it  * niston  radius/wavelength. 

AKOL  » k f.  ■ 2.0  * it  * rectangle  length/wavelength, 
o 

AKOW  =*  k _w  = 2.0  * it  * rectangle  width/wavelength. 

The  wavelength  is  of  F0  (F0  - (F1+F2) /2 . 0) . If  I0P  - 0 
and  I0D  = 0 then  this  card  must  be  omitted. 

SLS,  PFRI,  DB,  N*1 , ID 

SLS  ■ the  value  in  dB  of  SLO*  at  which  a pattern  is  desired. 
PFP'  ■ the  value  of  FO/F-  at  which  a pattern  is  desired. 


II 


l)H  ■ the  lower  limit  of  Ihe  pattern  plot  (e.g.  1111  -*  -80.0) 

Nf  = the  number  ol  tlu*t;i  value?;  Ml  wlileli  tin*  beam  pattern 
Is  evaluated. 

TD  ■ AO  *■  the  step  size  between  theta  values  in  degrees. 

If  10P  ■ 0 this  card  must  be  omitted. 

Some  of  the  sections  have  graphs  with  multiple  curves.  The 
program  presently  limits  the  number  of  curves  per  graph  to  five. 

The  graphs  that  have  SLO*  as  the  X-axis  are  presently  limited 

to  twentv-one  (21)  separate  values  of  SLO*.  These  are  generated  in  the 

% 

program  by  the  formula  SLO*  ■ 180.0  + 10.0  * (J-l),  J » 1,  2,  3,  . . . , 21. 

The  maximum  number  of  theta  values  for  patterns  Is  restricted  to 
101  because  of  the  size  of  the  plotting  grid. 

If  a pattern  Is  desired  from  a rectangular  source  the  pattern 
Is  assumed  to  he  in  the  plane  of  the  length  axis  of  the  rectangle. 

For  the  directivity  index  of  a piston  the  pattern  is  assumed 
to  be  zero  behind  the  piston. 

Three  subroutines  are  used  in  the  program — BMPAT,  DLPAT,  CNPAT. 
BMPAT  is  used  to  calculate  the  beam  pattern  for  IOC  * 1. 

DLPAT  is  used  to  calculate  the  beam  nattern  for  IOC  * 2. 

CNPAT  is  used  to  calculate  the  beam  pattern  for  IOC  = 3. 

The  program  calls  on  several  packaged  routines  to  evaluate 

one-dimensional  definite  integrals,  N'-dimensional  definite  integrals, 

"I"  Bessel  functions,  "J"  Bessel  functions,  and  exponential  integrals. 
Packaged  plotting  routines  are  used  as  well  as  a routine  to  determine 
the  amount  of  CPU  time  used  bv  the  program.  Ail  of  these  routines  are 
described  in  Anotndix  A. 
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I’ROCRAM  OUTPUT 

Figure  l shows  an  example  of  l lie  program  output  for  20  log(ETA), 

-5  2 

SrL-*,  and  SPL(OUT)  for  a range  of  u.,.ro  values  from  10  dB  to  10  dB. 

Five  values  of  FO/K-  were  used  varying  from  5.0  to  80.0. 

Figure  2 shows  an  example  of  the  program  output  for  the 
patterns  of  a piston  with  a kQa  of  10.0  where  FO/F-  - 10.0  and  SPL0*  Is 
varied  from  240  dB  to  330  dB. 

Figure  3 shows  the  program  output  for  the  directivity  Index  of 
a piston  with  a k^a  of  10.0  and  an  FO/F—  of  10.0. 
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APPENDIX  A 

Included  In  this  appendix  are  descriptions  of  the  packaged 
routines  that  are  used  in  the  program.  Tlu*  routines  ROMBS,  R0M2,  RMB1, 
RMB2,  RMB3,  BJY01,  BK'n,  and  1)1  I are  from  the  JPL  Fortran  V Subroutine 
Directorv,  Edition  No.  4,  October  1970.  The  routines  PSCA1.E,  PSETUP, 
PPLOT,  and  SECOND  are  from  the  Westinghouse  Research  Laboratories' 
Fortran  library. 
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11.1.  one-dimensional 

11.1.1.  QUADRATURE.  ONE-DIME US I ON,  S.P. 

11.1.1.1.  IDENTIFICATION 

QUADRATURE,  ONE-D I KENS I ONAL , SINGLE  PRECISION 

LANG  FILE  El.T/vERS  SI?E  ENTRY  NAmES 

F-V  LIB*JPL$  ROMPS/ JPL  1?S1  ( P.  > r.ftOl  ( IP  > ROMBS.ROM2 

SUBROUTINES  USED:  MNONEtt 

COGNIZANT  PERSONS:  W.  R.  BUMTON  AND  M.  OIETHFl.M' 

JPL,  SECTION  3l'i,  I960  SEPT  3f) 


11.1.1.2.  PURPOSE 

OBTAIN  APPROXIMATE  EVALUATION  OF  A ONE-DIMENSIONAL  DEFINITE 
INTEGRAL  UY  numerical  quadrature. 

ANS  =.  THE  INTEGRAL  FROM  A TO  D OF  F(XH=DX 


11.1.1.3.  P!  FERENCrS: 

FOR  A COMPLETE  DESCRIPTION  OF  THIS  SUBROUTINE , INCLUDING  A 
DISCUSSlOu  OF  A VARIETY  OF  TFST  CASES,  SEE: 

WILEY  R.  bUNTON,  MICHAEL  UlETHEl-M,  AND  KAREN  MAlGLER,  'ROMBERG 
QUADRATURE  SUBROUTINE  TOP  SINGLE  AND  MULTIPLE  INTEGRALS'.  JPL 

W,  E MCN,  M.  DIET"EL  '?  G.  WIUjE,  * MOD  I f TED  ROMtiFPG  QUADRATURE:  A 

SUBPCjTliJ-  TO  SUPPORT  GENERAL  SCIENTIFIC  COMPUTING*,  jPL  X uTEPflAL 
MS"CPTNGUM  T’  314-ESB-  A^PIL  1-  1 97j, 

W.  BUNTOIJ,  M.  DlETH"Lf- , 'MO”  IF  -r AT  IONS  To  THE  JPL  POMDFR& 

SUBROUTINES' f 7M  31;:-c47,  I i07i:. 


11.1.1.4.  METHOD 

THIS  SUBROUTINE  COMBINES  TECHNlO"ES  FROM  'ROMBERG'  AND  'AOAPTIVE 
STEP'  QUADRATURE  METHODS.  THE  SUBROUTINE  TNITIAlLY  PICKS  A 
SUBXNTERVAL  [A, BIT  OF  THE  TOTAL  TNVFRVAL  f A,H]  AND  attempts 
TO  APPROXIMATE  THE  INTEGRAL  OVER  T!  IS  S^'INTEP'/AL  PY  USING  THREE 
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STAGES  OF  ROKIIFRG  QUADRATURE.  U1IS  RioUIRf  S EVALUATION  OF  THE 

integrand  at  five  equally  spaced  points. 

Ir  THIS  QUADRATURE.  IS  REGARDED  AS  SUCCESSFUL  (SEE!  ERROR  CONTROL) 
THE  SUBROUTINE  WILL  ADO  TK!  VAl.UF.  OBTAINED  TO  A RUNNING  SUM  AND 
PROCEED  TO  TREAT  A NEW  DISJOINT  SufUNTL'RvAL  OF  THE  SAME  OR 
GREATER  LENGTH, 

IF  THIS  QUADRA 'URL  IS  |\|OT  REGARDED  AS  SUCCESSFUL  AND  THE  CURRENT 
STEP  LENGTH  IS  CRLATER  TUAN  HfUN*  THEN  THE  SUBROUTINE  WILL 
REJECT  THE  RESULT  FOR  THE  C UP RE NY  SUB INTERVAL  AND  TAKE  THE  LEFT 
HALF  OF  THE  SUblNTERVAL  AS  Till"  HEW  SUUlNTf RVAL  TO  RE  TREATED. 

IF  THE  CURRENT  STEP  LENGTH  IS  HMJM  OR  SMALLER*  THEN  THE 
SUBROUTINE  WILL  ACCEPT  THE  CURRl  UT  RESULT  AN!)  PROCEED  TO  THE  NEXT 
SUBINTERVAL»  WRlTJNG  A MESSAGE.  ON  FORTRAN  UNIT  (>  IDENTIFYING  TlC 
SUBINTERVAL  Oil  WHICH  VnE  ACCURACY  TFST  WAS  NOT  SATISFIED. 


11.1.1.5*  ERROR  CONTROL 

THE  QUADRATURE  IS  REGARDf D AS  SATISFACTORY  OVr P A PARTICULAR 
SUBINTERVAL  IF 

1,  THE  ESTIMATED  RELATIVE  ERROR  OE  THE  QUADRATURE  OVER  THAT 
SUB INTERVAL  IS  AT  MOST  E UMAX  r OP 

2 THE  ESTIMATED  ERROR  IN  THAT  Si  '0 1 NTFRVAL  RELATIVE  TO  THE 

ACCUMULATED  VAl’uE  OF  TllE  INTEGRAL  UP  TO  AND  INCLUDING  THAT 
SUBINTERVAL  IS  at  MOST  I'.RMAX. 

3.  THE  ESTIMATED  ABSOl  UTF.  ERROR  OVER  THAT  SUB  INTERVAL  IS  AT 
MOST  . 1 *AbS (ERmAX ) . 


11.1. I. 6.  PARAMETER  CHECKING 

the  subroutine  will  terminate  execution  with  a printed  message 

A JD  A STANDARD  EXEC  R * ALE -RAO'  (RETURN  C)  JE  ThL  GIVEN  PARAMETERS 
DO  NOT  SATISFY 


A<B»  Op  A>B » BUT  MOT  ArB 
0.<hmin.le,hstar.le.umax  AND 

RELATIVE  - 0.<ERMAX<1 ,C#  BUT  NOT  FPMAX=0. 
absolute  - any  negative  NUMBER*  BUT  NOT  ro. 

IF  A.GT.B  OR  IF  HSTAR.GE, (B-A)/4*  ONLY  A WARNING  MESSAGE 
IS  PRINTED. 
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11.1.1.7.  USAGE 

REAL  A t U.  X,  FOFX.  Mf,TAR»  M'llN,  UMAX » ERMAX » ANS 

integer  k » key 

[ASSIGN  VALUES  TO  A.  B,  Hf,TAR»  HMIN,  UMAX.  ERMAX,  AND  KpY  ] 

CALL  ROMtjS  ( A » D»  X .FOFX.MSTAR » HMIN .UMAX » ERMAX » ANS, K »KEY ) 

1?  [EVALUATE  THC  INTpGRAUO  USING  THE  CURRENT  VALUE  OF  X AND 
STORE  THt  RESULT  IN  T0FX.3 

CALL  RUK2 

IF(K.EQ.l)  GO  TO  10 

[QUADRATURE  IS  COMPLETED.  RESULT  IS  IN  ANSI 


THE  SUBROUTINE  PARAMETERS  ARE  DEFINED  AS  FOLLOWS: 

A.B  LIMITS  OF  INTEGRATION.  REQUIRE  A NOTrB. 

X VARIABLE  SET  BY  TMF  SUBROUTINE  FOR  INTEGRAND  EVALUATION 

IN  THE  USCR’S  PROGRAM. 

FOFx  VALUE  OF  INTEGRAND  COMPUTED  BY  USER’S  PROGRAM  USING  Tht 

ARGUMENT  X, 

HSTAR  SUGGESTED  INITIAL  STEP  SIZE.  The  INITIAL  STEP  SIZE,  H» 
WILL  BE  SET  AT  Hr.OK  <P~A>  IF  HSTAR. C»r  • (B-A)/4,  OR  AT 
HrHSTAR  IF  HSTAR.I  T.  (D- TF  F FIR‘.*i  SUBINTERVAL  WILl- 
BE  OF  LENGTH  tu*M  AND  WILL  REQUIRE  FIVE  EVALUATIONS  OF 
THE  INTEGRAND  AT  X*A,  A*II,  , . , A*“tH.  SUGGEST 
HSTAR. GE.  (B-A)/i4 . 

REQUIRE  HM In  .LE.  HSTAR  ,LF.  HMAX. 

HMIN  MINIMUM  ALLOWABLE  ‘‘TCP  SIZE . Rf.OUIPt  0.  < HMJN  .LE, 

HSTAR 

HMAX  MAXIMUM  ALLOWABLE  STEP  SIZE,  REQUIRE  HSTAR  .LE.  HMAX. 

ERMAX  TOLERANCE.  ON  RELATIVE  OF  ABSOLUTE  rRf'OP.  SFC  DISCUSSION 
Arovr.  UNDER  ’ERROR  CONTROL’.  REASON  mill  settings  FOR 
tR"AX  WOlll  D BE  IN  THE  EAnGE  FROM  lrF-<>  TO  l.E-7.  IF 

greater  ac curacy  ir  RrouiRfD-  sn  the  v rite-up  on  rombd. 
IT  is  REQUIRED  THAT  O.LT.FRV  ;,LT , 1.  FOR  THE  RELATIVE 
ERROR  TEST,  ONE  SHOULD  KNOW  TH£  RmMGE  OF  THE  ANSWER 
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ANS 

K 

KEY 


11.1 

1. 

2. 


BEFORE  me  uses  ABSOLUT » f RROH . 

THE  FINAL  VALUE  OF  THf  INTEGRAL . AVAILABLE  WHEN  ROM2 
RETURNS  WITH  K=2, 

BRANCHING  FLAG  SET  BY  THF  SUBROUTINE  FOR  USE  IN  THE 
USER'S  PROGRAM,  Krl  MEANS  ThE  USER  SHOULD  EVALUATE  THE 
INTEGRAND  AT  X,  STORE  THE  VAL"E  IN  FOFX.  AND  RE-ENTER 
R0M2,  Kr2  MEANS  THE  COMPUTATION  IS  COMPLETED  AND  THE 
VALUE  IS  IN  ANS. 

FLAG  TO  CONTROL  PRINTING  OF  ERROR  MESSAGES,  PREVIOUSLY , 
ANY  VALUE  OF  KFY  NOT  = 7 WOULD  WRITE  A DIAGNOSTIC  MESSAGE 
WHEN  H BECAME. LT.HMIN,  THE  INPUT  VALUES  HAVE  BEEN 
CHANGED  SO  THAT  WHEN 

ACTION 

KEY='j  PRINT  INTERNED  I ATI  T AND  Y VALUES  I PRINT  THE  HMlN 
DIAGNOSTIC  IF  DETFCTED. 

=G  PRINT  INTERMEDIATE  T AND  Y VALUES*  DO  NOT  PRINT 
THE  HMlN  DIAGNOSTIC. 

-7  DO  NOT  PRINT  TltE  T AND  Y VALUES  OR  HMlN 
DIAGNOSTIC. 


=A|jY  OTHER  VALUE  PRINT  THE  HMlN  DIAGNOSTIC  IF 

detected. 

THE  T VALUES  PRINTED  t* E T(1.CJ»  T(l,l)»  AND  T(2.0).  THE 
PRINTED  Y VALUES  APE  THE  FUNCTIONAL  EVALUATIONS  Y(l)  THRU 
Y ( 5 ) AT  THE  POINTS  X»X+H, . , ..X+4H,  SEE  REFERENCES. 


.1.8.  remarks: 

the  following  notes  may  be  applicable  for  difficult  integrands. 

IF  IN  DOUBT,  One  SHOULD  USE  A SMALL  VALUE  nr  HSTAR.  THE  STEP 
SIZE  H CAN  DOUBLE  QUICKLY  AND  THE  USER  IS  PENALIZED  ONLY  A 
SMALL  number  On  FUNCTIONAL  [VALUATIONS  WHILE  HE  INCREASES  HIS 
CHANCES  OF  GETTING  AN  ACCURATE  APPROXIMATION  MANY  FOLD. 

BE  CAUTIONS  WHEN  RELATIVE  ERMAX  IS  .GT.  in-j.  IF 
HSTAR.LT.  Ui-A)/L’»  BUT  NOT  SMALL  ENOUGH,  «ND  THF.  FUNCTION  IS 
OCC1LATORY,  VERY  DIFFICULT,  ETC.,  HOMDS  CAN  RETURN  A WRONG 
ANSWER,  EXAMPLE!  F(V»=X*SlN3nXtCOSX  ON  THE  INTERVAL  (0.2PI), 
HSTaR-1 . 57 , TRUE  ANSWC 20^672 40.  A RELATIVE.  TOLERANCE  OF'  .1 
GAVE  -.25L-S*  AND  A RELATIVE  TOLERANCE  OF  .{'I  GAVE  4.1801  BAO 
ANSWERS. 
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3,  ALSO  FOR  RELATIVE  FRMAX  ,LT,1C**-G.  IF  THE  Rr LATIVE  ERROR  I S 
ASKIMG  FAR  GREATER  THAN  (•>  SIGNIFICANT  DIGITS  ONI  IS  PUSHING  THE 
ACCURACY  <>F  THE  UNI VAf  11CH.  ON  A HIGHLY  OCHiATOIIY»  VERY 
DIFFICULT  PHOHLEM.  ROMMS  hay  hi:  TAKING  thousands  morl 
FUNCTIONAL  EVALUATIONS  and  Ijor  RrAUdNG  the  accuracy  it  did  at 

example:  same  f<xi  as  amove  when  gjvi  n tourancf  wa  > 

10*0-6#  AMS=..7Cc>673nO  WITH  A .’OCO  HlfICTlOilM  I VAIIAT IOllS » HUT 
V.HCN  WHEN  ERMAX=10**-n»  ANS  •«'0*l>776f»  AND  10 » 267  FUNCTIONAL 
EVALUATIONS. 

4,  IF  ONE  WANTS  A ROUGH  APPROXIMATION  OF  THE  IflTRGRAL#  HE  CAN  SET 
HMIN=HSTAK=HMAX . A Fixrn  ritrP  IHTRGRATION  01-  THE  FUNCTION  WILL 
TAKE  PLACE.  BL  SUPf.  THAT  I'l.'-V*  OP  MANY  [)l AbNOST ICS  MAY  BE 

printed. 
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11. r.  H»»LTI-r»IMENS10rjA! 


11.2.1.  OUADUATURl#  MUl.TI-niMENf.lOM,  S.P. 


11.2.1.1.  inf  NT  IF I CAT I ON 
QUADRATURE,  MULTI-DIMENSIONAL,  SINGLE  PRECISION 


ENTRY  NAME 
RMB 1 

RMBlA,f<MH2 
RMB  3 


SUBROUTINES  USED:  MNONEW 

COGNIZANT  PERSONS:  W.  R.  Rt INTON  AMD  M.  DlETHELM. 

JPL , SECTION  31*1  • 1969  SEPT  3q 


LANG 

F-V 

F-V 


FILE 

LlbtJPLS 

LIb*JPL«- 


ELT/vEPS 
RMB 1 /JPL 
PMfi/JPL 


SI?L 

361 in»"24l ( 1°) 
2/67 ( h } - 1 207 (1C) 


- CALLING  sequence  to  this  subroutine 

- HAS  BEEN  CHANGED  SINCE  ED,  3 1 APRIL  70- 

- VERSION  OF  THIS  DIRFCTOPY. 


11.2.1.2.  PURPOSE 

OBTAIN  APPROXIMATE  EVALUATION  OF  AM  N-plMFNSlONAL  DEFINITE 
INTEGRAL  UY  NUMERICAL  QUADRATURE.  THE  LIMITS  Of  INNER  INTEGRALS 

can  be  functions  of  thf  variables  of  the  outer  integrals. 


result  = THE  INTEGRAL  FROM  A1  TO  B1 
with  respect  to  y < i > or 

ithe  integral  from  n to  nr  with  pfe^f r t to  xc2)  of 

C...CThE  INTEGRAL  FROM  AN  VO  nr-  V/ J * f'FSPFC.T  TO  X(N)  OF  Fl 

...33 

WHERE 

1.  Ml  AND  PI  ARC  CONSTANTS 

2.  FOR  J=2»...#N-  AJ  AND  BJ  uAy  PF.  FUNCTIONS 
OF  X ( 1 ),..., Y ( J-  ) , 


3.  THE  INTEGRAND  F, 


IS  A FUNCTION  OF  X ( 1 > » . . . , X (N) . 
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11.2.1.3.  RtFE.RF.NtFC 

FOR  A COMPLETE  DESCRIPTION  OK  THIS  M(,iPOuT  iNf. , INCLUDING  A 
DISCUSSION  OF  A VAUIMy  Or  1 1*  ST  OASIS,  SF  E : 

WILEY  k.  UUNTONt  MICHAEL.  DIF  TUf  l.M,  ANN  K A;  T I I HAK.LFR,  • ROMflERG 
QUADRATURE  suijroutine  for  single  and  multipif  integrals*  ♦ JPL 
INTERNAL  MEMORANDUM  TM  3l»4'*22l»  JULY  If  Ifl6'J.  '-Ft  OTHER 
REFERENCES  UNDER  ROMOS. 


11.2.1.0.  METHOD 

THE  SUBROUTINE  PERFORMS  A NESTED  SR’ULNCt  Of  ONE  ~0  IMF.  NSIONAL 

integrations,  thl  suuroutime.  usrr  iti:  the  oni -dimensional 

IllTEGK.  TIOIIS  IS  ESSENTIALLY  ROM'BS  U/TU  CIlMNGf  S IN  THr  WAy 
STORAGE  IS  managed.  stf  the  WRlTf  -UP  ON  I’OMB'j . 


II. 2. 1.5.  USAGE 

THE  USER  MUST  PROGRAM  All  H II T I Al  I7ATI  ON  CALL  TO  Rvni,  THEN  N 
SEPARATE  L-'LLS  TO  RMB?,  ONE  for  r/Cil  LEVEL  of  INTEGRATION.  AND 
FINALLY  A CALL  To  RMBJ , SOMC  01  TtifSF  CALL  STATEMENTS  W 1 1 L DE 
executed  MORE  than  onci  . 

THE  CODING  SPECIFICATIONS  ARf  AS  FOLLOWS: 

integer  ii.kgo.kfy 

real  x i iii i 1 1! ) t r,  result,  fpmax*  woiuan) 

REAL  Al.  Ul » IlSTAIUf  iPUfU.  HIlAxi 


REAL  an,  dm,  hSTAPN,  hMInN,  mmaxn 
warning:  calling  scquencf  to  tuts  suproutine  has  been  changed 

SINCE  ED.  :.  1 APRIL  lovo  vmSJON  cr  THIS  DIRECTORY. 


I ASSIGN  VALUES  TO  N AND  TRmAx.1 
CALL  RMD1  (M,X.F  ,f  [ r I *L  T , ERpAx  , KGO , W » K[  Y) 

1 l ASSIGN  VALUES  to  all  PARAMETERS  tn  the  FOLLOWING 

CAI  L STATEMENT.] 

CALL  RKO£  ( A I r n 1 .IlSTARl , HM T ’ U , I IVAX  1 ) 

2 t ASSIGN  VALUES  TO  ALL  PARAMETERS  IN  THL  r Ol.LOWlNG 

CALL  STATEMENT,  SOmF  O'  Tl  F?E  VALUES  MAY  BE  COMPUTED 
AS  FUNCTIONS  OF  X ( I } * i 
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r ASSIGN  VALUES  TO  ALL  RARAMF.TF.RS  IN  THE  FOLLOWING 
?ALL  STATMlrNT.  SOME  OF  THE  VALUES  MAY  BE  COMPUTED  AS 

FUNCTIONS  OF  X(  i )»..*■  »X(M-1  )•  J 
CALL  RRBL  ( AN . nri  r MSTARN f HmINN  f I IMAXH) 

Riain  L COMPUTE  F as  a function  of  

CALL  RMH: 

GO  TO  (l»2»...»WN«#l»N*lM,«N*rrt?#K60 

LTHE  COMPUTATION  IS  COMPLETED.  THE  VALUE  OF  THE  INTEGRAL 
IS  IN  RESULT. 1 


MN+PH 


THE  DlMENSIOi  1 1 f-jG  PARAMETERS  MUST  SATISFY; 

Ill'll « .OF..  N 

«n?m  ,gf.,  ron: 


THE  PARAMETERS  FOR  PMR1  APF  DEFINED  AS  FOLLOWS; 
N 


NUMF:ER  OF  LF  VrLS  OF  INTEGRATION  REQUIRED.  (N  .GE.  1) 
EXAMPLE;  For  a TRIPLE  INTEGRAL*  SF.T  NrB. 


(Kill  f'lti'll  CURRENT  VA1  urf  OF  TuF  INTEGRATION  VARIABLES  SET  RY 
THE  JlHROUTIIIl.  X<  1 > IS  ASf.OClAIED  WITH  THE  OUTERMOST 
integration  and  y < m ) witu  the  t-nfrmost. 

F VALUE  OF  INTEgRAIH).  to  nr  commuted  by  the  user  as  a 

FUNCTION  Or  x ( 1 ) » . • . > * ( U ' AFTf  R TFT  N-TFF  CALL  T 
RMnP  AND  WHEN  Rflfi  RETURNS  WITH  I'.GO  = N+l. 

RESULT  VALUE  OF  I NTT  ORAL  COMPUTET  RY  THF  SUBROUTINE.  AVAILABLE 
WHEh  RMlib  RETURNS  VT.Tr*  KGO  - N*P. 

ERMAX  RELATIVE  CR  ABSOLUTE  ERROR  TRLEtaNCE  SET  RY  THE  USER. 

RFASONABi  F VA>  JtS  WOU1  n RE  a'i  THE  RANGE  l.E-P  To  l.E-4 
FOR  RLlATIVL.  REQUEST ING  '*ORF  ACCURACY  MAY  GREATLY 
INCREASE  THE  EXECUTION  TIME.  FPMAX  MUST  SATISFY 

C .lt.ermax.lt , i , ron  relative*  anv  negative  number  not  = 

c FOR  ABSOLUTE . 

KGO  branching  flag  sft  nY  the  subroutine  rmb?  for  use  in 

THE  USER’S  PROGRAM. 

KGO  r I,..., IJ  MEANS  THE  USER  SHOULD  SE-T  THE  PARAMETERS 
for  The  KGO-TH  call  TO  rmbp  and  then  EXECUTE  THAT 


\ 
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CALL. 

KGO  = N+l  MEANS  THE  USER  SHOULD  COMPUTE  THE  INTEGRAND* 

F#  AND  THEN  CALL  RM113. 

KGO  = N+2  MEANS  THE  COMPUTATION  If,  COMPLETED  AND  THE 
VALUE  OF  THE  INTEGRAL  IS  IN  RESULT. 

(MCI). I=1 »29*M)  WORKING  SPACE  NEEDED  OY  THE  SUBROUTINE. 

KEY  FLAG  TO  CONTROL  PRINTING  OE  ERROR  MESSAGE  WHEN  ERROR 

TOLERANCE  IS  NOT  MFT  WITH  MINIMUM  STEP  SIZE. 

KEY=7  NO  ERROR  MESSAGE  IS  PRINTED 
KEY=AhY  other  VALUE  error  message  is  printed. 

THE  PARAMETERS  FOR  THE  J-TH  CALL  TO  RM02  ARE  AS  FOLLOWS! 

AJ»BJ  LOWER  AND  UPPER  LIMITS  OF  INTEGRATION  FOR  THE  J-TH 

INTEGRAL.  E TTHER  AJ  > DJ  OR  AJ  .LE.  BJ  IS  PERMITTED. 

HSTARJ.HMINJfHMAXJ  QUADRATURE  STFP  PARAMETERS  FOR  THE  J-TH 

integral,  these  PARAMETERS  MUST  BL  POSITIVE*  REGARDLESS 
OF  THE  SIGH  OF  <f!J  AJ)r  AND  MUST  SATISFY  0.  < HMINJ 
,LE . HSTARJ  .LE.  HMAXJ.  THE  J-TH  QUADRATURE  WILL  START 
*»ITH  A STEP  WHOSE  MAGNITUDE  IS 

HJSMINCHSTARJ, AHS(BJ-aJ)/4. ) EXCEPT  FOR  OUTER  INTEGRAL. 
SEE  REMARK  4.  BELOW. 

AND  THE  MAGNITUDE  OF  THE  STEP  WILL  BE  KEPT  BETWEEN  HMINJ 
AND  HMAXJ. 


11.2.1.6.  REMARKS 

1.  THC  VALUE  OF  KGO  WILL  NEVER  BE  1.  THUS  THE  FIRST  CALL  TO 
RMB2  WILL  ONLY  BE  EXECUTFO  ONCE, 

2.  THE  PARAMETERS  N»  ERMAX.  Al'D  AI  L PARAMETERS  IN  THE  CALLS  TO 
RMB2  WHICH  ARE  70  REMAIN  CGHETANT  DURING  THE  QUADRATURE  CAM 
BE  SPECIFIED  l-ITEFA|  I.Y  IN  7!  CALL  STATEMENTS, 

3.  PARAMETERS  IN  THE  SECOND  THROUGH  THE  N-TM  CALL  TO  RMR? 

WHICH  ARE  NOT  CONSTANT  CAN  pC  WRITTEN  IN  THE  CAI  L STATEMENTS 
AS  ARITHMETIC  EXPRESSIONS.  THE  PARAMETER  E*  HOWEVER  * 

CANNOT  BE  REPLACED  i N THE  CALI.  TO  RMRl  BY  An  ARITHMETIC 
EXPRESSION  EXCL.PT  FpR  THE  SPECIAL  CAFF  In  WHICH  r Ic 
CONSTANT. 

4.  IF  THE  STARTING  STEP  SIZE  FOR  THE  OUTER  INTEGRAL  IS  GREATER 
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THAN  OR  FOIIAL  TO  < M-A  )/«i.*  A WAI.K-RAFk  DIAGNOSTIC  (RETURN  0)  IS 
PRINTED  AND  rXlf.llTIOII  IS  TERM? NATEC. 

5,  SUGGESTED  VALUE  FOR  Kl. Y ID  G, 


11.2.1.7.  EXAMPLE 

CASE  11*  PAGE  4«*  OF  THE  REFERENCE  CITED  ABOVE  CAM  BE  CODED 

AS  follows  : 


Rt AL  X( 3) * WORK ( R7 ) 


CALL 

RMS  1 ( 3 , X , F , RFSUL T , ! . E-4 , KGO , 

WORK , 0 ) 

1? 

CALL 

RMU£ ( P , , 

t . , 

l.F-2,  l.E-4, 

1.) 

20 

call 

rmoe ( X ( 1)  * 

X ( 1 ) + * 2 , 

i 0 E-2  * 1.1  —4  , 

1. ) 

S'1 

CALL 

Ri-Ujt  ( X ( ] )*X(2)  , 

X ( 1 ) +X ( 2 ) , 

l.E-2,  1 .1-4, 

1.) 

4 0 F r x(l)*X(2  XX(o) 

CALL  RMbl 

GO  TO  (10,  2C  * 30,  ur,  SO)*  KGO 
S?  CONTINUE 


THE  TRUE  VALUE  OF  THIS  IllTTGRAl.  IS  -C . C 3.200(1 1 OS  . THE  COMPUTED 
VALUF  HAD  AtJ  ABSOLUTE  ERROR  OF  S.1F-?  A'Hj  A Rf  L AT  I VC  ERROR  OF 
1.6E-7  . TM  S"BRO(JTIIlF  WROTE  Two  DIAGNOSTIC  MESSAGES  INDICATING 
THAT  ThE  REQUESTED  ACCURACY  COlil.t)  NOT  Elf  ATTAi:,..D  WITH  THE 
MINIMUM  STEP  SI2E  0!  i TnF  INTERVALS  f . C J.F  . x*l)  .LE.  4.E-4  AND 
4.E-E  .LE . X ( 1 ) ,LL.  8.r--«!  , Tills  EXAMPLE  USED  10361  EVALUATIONS 
OF  TUB  IIITLGivAijD, 

THIS  EXAMPLE  WAS  RUN  2c  TIMES  AND  TIMFD  USING  PRTIM1/PRTIM2. 

THE  TIME:  VARIED  from  2. Gil  SECONDS  TO  r.,943  seconds. 
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‘*.2.6.  BESSEL  JC.J1.Y0.Y1.  O.P. 

‘♦.2.6.1.  IDEi'IT  IFICATIOKI 

bjy?i/double  precision  bessel  functions  jo.  ji*  yo»  and  yi 


LANG  FILE  F.LT/VERS 

F-V  LIB*JPL*  flJYOl/JPL 


SIZE  ENTRY  NAME 

iniu(8)=52lM10)  BJY01 


SUBROUTINES  USED:  DSQRT,  OSIN,  DCOS.  OLOG.  AND  CHBPOL 

\ 

COGNIZANT  PERSON!  E.  W.  NG.  JPL.  SECTION  314,  1969  AUGUST  1 
*♦.2. 6. 2.  PURPOSE 

TO  COMPUTE  THE  VALUES  OF  THE  BESSEL  FUNCTIONS  JO,  Jl,  YO#  ANo  Yi. 
4.2.C..3.  ACCURACY 

FOR  X .LE.  N,  AT  LEAST  15  SIGNIFICANT  DIGITS* 

FOR  X > N»  AT  LEAST  IS  OECIMAL  PLACES. 


4.2.0. 4,  USAGE 

Integer  no.  ru 

DOUBLE  PRECISION  X,  RC<2>»  B1C2) 

CALL  BJYOKX.Nn.Nl.BD.Bl) 

THE  PARAMETERS  are  defined  AS  follows: 

X ARGUMENT  AT  WHICH  FUNCTION  EVALUATION  IS  DESIRED. 

REQUIRE  x > n.  IF  X .LE.  U THE  SUBROUTINE  WILL  PRINT 
ON  FORTRAN  UNIT  l THE  MESSAGE: 

FOR  X=0.»  J0=1.  J1=0.»  YC=Y1=-INFINITY 

AND  WILL  SET  B0<1)=1.’  B1(1)=0.,  B0(2)sBl(2)=-l.D3B. 


NO 


n: 


=0  00  NOT  COMPUTE  jC  OR  YO 

=1  COMPUTE  JO  BUT  NOT  YU 

s 2 COMPUTE  JO  ANO  Yn 

=0  00  NOT  COMPUTE  Jl  OR  Yi 
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=1  COMPUTE  J1 

BUT 

NOT 

=2  COMPUTE  J1 

AND 

Y 1 

B0(  1 ) 

COMPUTED  VALUE  OF 

JO 

BO  ( 2 ) 

COMPUTED  VALUE  OF 

YO 

Bl(l) 

COMPUTED  VALUE  OF 

J1 

BI  (2) 

COMPUTED  VALUE  OF 

Y1 

4. 2. 6. j.  REMARKS 


THIS  SUBROUTINE  WAS  ORIGINALLY  DESIGNED  FOR  POSITIVE  X.  IF  IT  IS 
DESIRED  TO  INCLUDE  ITS  UTILITY  FOR  NEGATIVE  X AS  WELL*  THE 

fu.lowing  identities  should  be  used: 


Jp(-X)  = J0(X) * Jl(-X)  = -J1 (X) 
Y?(-x)  = Y0(X)  ♦ 2*SQRT < -1 ) * JC ( X ) 
YK-X)  = -Y1  (X)  ♦ 2#S0RT(-1UJI(X» 


note  that  y of  negative  argument  is  in  general  complex. 

(REFERENCE!  NbS  HANDBOOK  OF  MATH»L  FUNCTIONS*  APPLIED  MATH  SERIES 
NO.  55*  1964,  P.  36l> 
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4.2.9.  BESSEL  I AnO  K OF  GENERAL  OROER,  O.P. 

4. 2. 9.1.  IDENTIFICATION 

DOUBLE  PRECISION  BESSEL  FUNCTIONS  I AND  K OF  GENERAL  ORDER  NU 
AND  ARGUMENT  X. 

LANG  FILE  ELT/vERS  SIZE  ENTRY  NAMES 

F-V  LIB*JPL»  OESI/JPL  17?4(8)=1020(1C>  r<ESI»  BESK 

SUBROUTINES  USED:  OSIN#  DLOG#  OEXP 

COGNIZANT  PERSON!  E.  W.  NG.  JPL,  SECTION  314,  1969  JULY  23 

4. 2. 9. 2.  PURPOSE 

TO  COMPUTE  IN  DOUBLE  PRECISION  ARITHMETIC  THE  VALUES  OF  BESSEL 
FUNCTIONS  I AND  K OF  THE  GENERAL  ORDER  NU  AND  ARGUMENT  X. 


4. 2. 9. 3.  METHOD 

RECURRENCE  TECHNIQUES  AND  THE  PHASE -AMPLITUDE  METHOD  ARE  USEO. 
(SEE  REFERENCES! 


4. 2. 9. 4.  ACCURACY 

FOR  NU  .GE,  X,  A MINIMUM  ACCURACY  OF  lb  SIGNIFICANT  FIGURES  IS 
FOUND  I FOR  NU  < X,  A MINIMUM  ACCURACY  OF  15  DECIMAL  PLACES  IS 
INSUREO. 


4. 2. 9. 5.  USAGE 

OOUBLE  PRECISION  X,  V,  A(LDIMA),  AK(LDlMA) 

INTEGER  N,  LOIMA 

CALL  BESKX'VtN.A'LOIMA) 

OR  CALL  BESK ( X , V » N # A » AK , LD I MA I 

THE  ENTRY  BCSI  COMPUTES  THE  FUNCTION  DEXP(-X»  * I OF  ARGUMENT 
X AND  OF  ORDERS  V,  V+1,...,V*N. 

THE  ENTRY  BESK  COMPUTES  ROTH  THE  FUNCTIONS  DEXP(-X»  • I AND 
OEXP(X)  A K OF  ARGUMENT  X ANO  OF  THE  ORDERS  V»  V*1,...,V*N. 
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THE  SUBROUTINE  PARAMETERS  ARE  OEFJNEO  AS  FOLLOWS: 

X THE  ARGUMENT  AT  WHICH  FUNCTION  EVALUATION  IS  OESlREOl 

X .GE.  0.00 

V THE  NON- INTEGER  PART  OF  THE  OROER  OF  THE  BESSEL 

FUNCTIONS?  0.00  .LE.  V .LE.  1.00 

N IM4V  IS  THE  HIGHEST  ORDER  OF  TH£  BESSEL  FUNCTIONS 

OESlREOl  N MUST  BE  A POSITIVE  INTEGER. 

A(LOIma)  AIK).  K=l, , , . , (N4l)  ARE  The  LOCATIONS  in  which  the 

VALUtS  FOR  THE  BESSEL  FUNCTION  DEXP(-X)  * I OF  ORDER 
V#  (V*l),...,  <V«N)  ARE  RETURNED.  A ( ) MUST  HAVE 
DIMENSION  ,GE.  LOIMA 

AK(LOImA)  AMK),  Kri,...,  (N+i)  ARE  THE  LOCATIONS  IN  WHICH  THE 
VALUES  FOR  THE  BESSEL  FUNCTION  DEXP(X)  * K OF  THE 
ORDERS  V.  (V+l)#.,.#  (V^N)  ARE  RETURNED.  AM  ) MUST 
HAVE  DIMENSION  .GE.  LDlMA 

LOIMA  DIMENSION  OF  A ANO  AKl  IT  MUST  BE  *T  LEaST 

MAX(lOlNTCX)fN)  ♦ ZC  OR  MAX( lOINT (2**) #N)  FOR  X > 20, 


<1,2. 9. 6.  REMARKS 

FOR  X,  V,  OR  N BEYOND  ITS  REQUIRED  RANGE  OR  LDlMA  NOT  LARGE 
ENOUGH,  THE  MESSAGE  ‘ERROR  IN  BESSEL*  ANO  THE  VALUES  OF  X,  V» 
N#  NU,  ANO  LOIMA  ARE  PRINTEO,  THE  NUMBER  <N*2)  IS  THE 
DIMENSION  OF  a NEEDED.  IN  THE  CASE  X OR  V IS  OUT  OF  RANGE. 
THE  PARAMETER  NU  IS  irrelevant. 


4. 2. 9. 7,  REFERENCES 

M.  GOLDSTLIN  ANO  R.M.  THALER,  ‘RECURRENCE  TECHNIQUES  FOR  THE  CAL- 
CULATION OF  BESSEL  FUNCTIONS*#  MTaC#  VOL,  XIII,  PP.  102-108. 

H«,S?V°STEIN  AN°  f<*  M*  THALER»  ’BESSEL  FUNCTIONS  FOR  LARGE  ARGU- 
MENTS*, MTAC#  VOL.  XII,  1958#  PP.  18-26. 

• 

D.  JORDAN,  ARGONNE  NATIONAL  LAB*#  PROGRAM  WRITE-UP  NO,  C3715#  NOV 
1967.  * 
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4.2.13.  EXPONENTIAL  INTEGRAL 


4.2.13.1.  IDENTIFICATION 


exponential  integral 

LANG  FILE  ELT/VERS  SIZE  ENTRY  NAMES 

F-V  LIB* JPL*  DEI /JPL  DEI 


SUBROUTINES  USED:  DEXP.DLOG 


COGNIZANT  PERSONS  E.W.  NG,  JPLi  SECTION  314,  1970 » SEPT.  14 


4.2.13.2.  PURPOSE 

COMPUTE  THE  EXPONENTIAL  INTEGRAL  IN  DOUBLE  PRECISION  ARITHMETIC. 
FOR  X .GT.  0*  THE  EXPONENTIAL  INTEGRAL,  El,  IS  DEFINED  AS 

EI(X>  = INTEGRAL  FROM  Tr-INFINITY  TO  T=X  OF 
(EXP(T)/T»*0T 

WHERE  THE  INTEGRAL  IS  TO  BE  INTERPRETED  AS  THE  CAUCHY  PRINCIPAL 
VALUE.  FOR  X .LT.  0,  El  I XI  =-El(-X),  WHERE 

EKZIsINTEGRAL  FROM  T=Z  TO  t=infinity  OF 
<EXP(-Tl/T|*DT 


4.2.13.3.  METHOD 

THIS  SUBROUTINE  COMPUTES  THF  EXPONENTIAL  INTEGRAL  BY  CHEBVSHEV 
RATIONAL  APPROXIMATIONS  FROM  W.J.  CODY  AND  H.C.  THACHERf  JR., 
MATH,  COMP,  VOL,  22,  PP,  G41-65C,  AND  VOL.  23,  PP.  289-303. 


4.2.13.4.  ACCURACY 

EXTENSIVE  TESY',  WERE  PERFORMED  ON  THE  UNIVAC  1108  AND  THE 
FOLLOWING  accuracy  statistics  WERE  FOUNO: 

INTERVAL  MAXIMUM  RMS 

OF  X RELATIVE  RELATIVE 

<-150»  -4)  2.20-16  5,10-17 


r 

I 

!J 
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; U 

(-4#  -1)  7.10-17  1.20-17 

<-i»  0)  .4,70-18  1.40-18 

<?#  0.5)  1.80-16*  3 • L 0-17* 

(?.ii # 6)  5.50-17  1. CO-17 

(6»  12)  1.60-17  3.00-18 

Hit  24)  2.90-17  7.10-18 

(24.  IOC)  8.90-17  1.90-17 

CF.  E .to*  NG#  COMM.  ACM  VOL.  13»  «7.  PP.  448-449. 


4.2.13.5.  USAGE 

double  precision  OEI,  X 

USE  THE  FOLLUtolNG  FUNCTION  IN  A FORTRAN  ARITHMETIC  STATEMENT? 
DEI (X) 

WHERE  UEI  IS  THE  COMPUTED  VALUE  OF  THE  EXPONENTIAL  INTEGRAL. 


4.2.13.6.  LIMITATIONS 

CM0)=-IHFINITY  IS  APPROXIMATED  BY  -7.2D75  AMD  EI(X  .GT. 
174.673)  IS  APPROXIMATED  BY  *7.2D75. 

THESE  LIMITATIONS  WERE  ORIGINALLY  IMPOSED  FOR  THE  IBM/360  AND  HAVE 
NOT  BEEN  MODIFIED  FOR  THE  UNIVAC  iiu8. 


4,2.13.7.  ERROR  EXITS 
NONE 
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PPLOT 


Reference; 


D.  P.  Wei,  "A  Printer  Plot  Package  for  the  Univac  U1106", 
Research  Memo  72-1K4-COMPS-M16  (9-8-72). 


Purpose: 

This  package  of  subroutines  is  used  to  generate  graphs  on  the 
line  printer.  It  produces  an  8 by  10  grid  plot;  a grid  is  6 print 
characters  high  and  10  characters  wide.  Data  points  are  plotted  as  a symbol 
(e.g.,  *,  A,  3)  which  is  specified  by  the  user.  One  or  more  plots  can  be 
included  in  a single  graph. 

Usage: 

The  package  consists  of  three  Fortran  V subroutines. 

1.  PSCALE  scans  data  for  the  minimum  and  maximum  values, 

and  establishes  a range  of  values  which  produces  a rational 
scale  for  the  axis.  The  scale  is  chosen  such  that  each  grid 
is  1,  2,  or  3 times  an  appropriate  power  of  10. 

2.  PSETUP  sets  up  the  work  area  for  the  plot.  The  scales 

for  the  abscissa  and  ordinate,  and  the  labelling  information 
. for  the  axes  are  stored.  PSETUP  must  be  called  to  initialize 
each  graph. 

3.  PPLOT  stores  plot  data  and  generates  the  plot.  An  option 
controls  initiation  of  plotting.  Thus,  multiple  calls  on 
PPLOT  can  be  made  to  store  dat3  for  several  plots  before  a 
request  for  plotting  is  made. 


The  parameter  specifications  for  each  subroutine  follow. 


I 
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PSCALE 


Suhrout ine  Sncci f icst ion • 

SUBROUTINE  PSCALE  (V,  NPTS,  VLOW,  VHICH,  FIRST,  IXORY) 


DIMENSION  V(l) 
INTEGER  NPTS,  IXORY 
LOGICAL  FIRST 
REAL  VLOW,  VHIGH 


E * To  scan  the  data  in  vector  V , and  obtain  the  minimum  and 

maximum  values.  These  values  are  then  used  to  establish  a range  whir 
produces  a rational  scale  for  the  plot  axis. 


Usage: 


The  Fortran  calling  sequence  is: 

CALL  PSCALE(V , NPTS,  VLOW,  VHIGH,  FIRST,  IXORY) 


where: 


V ■ REAL  vector  containing  data  to  be  scanned. 

NPTS  ■ number  of  points  in  V to  be  scanned. 

VLOW, VHIGH  * output  variables  into  which  the  routine  returns  the 
lower  and  upper  values  for  the  range  of  V. 

FIRST  * LOGICAL  variable  or  value  indicating  when  it  is 
• the  first  set  of  values  to  be  scanned.  The  value 

is  true  (.TRUE.)  if  there  is  a single  vector  to  be 
scanned.  For  multiple  vectors  (when  more  than  one 
plot  is  to  be  included  in  a graph),  a call  on 
PSCALE  mu.it  be  made  for  each  plot  (vector).  The 
first  call  will  have  the  value  .TRUE.;  for  sub- 
sequent calls,  the  value  will  be  .FALSE.  The  range 
defined  bv  the  final  values  of  VLOW  and  VHIGH  will 
include  all  the  values  scanned. 

IXORY  - indicator  whether  abscissa  ra^pe  (x)  or  ordinate 

range  (v)  is  to  be  established.  Specifv  the  value 
0 for  abscissa  and  1 for  ordinate. 


' The  abscissa  axis  is  10  grid  wide.  The  ordinate  axis  is  8 gi id 

high.  PSCALE  need  not  be  used  if  the  range  of  values  for  the  abscissa  (or 
ordinate)  is  known.  The  number  of  grid  positions  should  be  taken  into 
consideration  for  proper  annotation  c f the  axis. 


1 

\ 

' 


I 
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PSETUP 

Subroutine  Specification: 

SUBROUTINE  (XLOW,  XHIGH , YLOW,  YHIGH,  IEOF,  ICRID,  LBLX, 

LBLY) 

REAL  XLOW,  XHIGH,  YLOW,  YHIGH 

INTEGER  IEOF,  ICRID 

DIMENSION  LBLX(l) , LBLY (1) 

Purpose: 

To  setup  work  area  before  data  values  to  be  plotted  are 
stored  (via  PPLOT  routine).  Labelling  information,  ranges  of  values 
for  the  abscissa  and  ordinate  variables,  and  options  for  grid  and 
annotation  are  supplied  as  input  parnu  eters . 

Usage : 

The  Fortran  calling  sequence  is: 

CALL  PSETUP  (XLOW,  XHIGH,  YLOW,  YHIGH,  IEOF,  ICRID, 

LBLX,  LBLY) 

where: 

XLOW, XHIGH  - the  lower  and  upper  values  for  the  range  of  the 
abscissa  variable  (x) . 

YLOW, YHIGH  * the  lower  and  upper  values  for  the  range  of  the 
ordinate  variable  (y) . 

IEOF  - Hollerith  option  to  set  either  E or  F format  for 

annotating  the  axes.  Specify  the  Hollerith  constant 
1HE  for  E-format  (1PE10.3),  or  1HF  for  F-format 
(FLO. 3). 

ICRID  - Hollerith  option  to  select  grid.  Specify  the 

Hollerith  constant  1HG  for  grid;  otherwise,  specify 
lHb  (b  denotes  a blank  or  space). 

LBLX  » INTEGER  vector  containing  information  for  label- 
ling the  x-axis.  The  first  element  LBLX(l)  contains 
the  number  of  characters  in  the  label.  The  labelling 
information  is  stored,  six  characters  per  word, 
starting  in  the  second  element  LBLX(2).  The  maxfit.T 
label  length  is  30  characters. 

LBLY  * INTEGER  vector  containing  information  for  labelling 
the  v-axis.  The  rules  for  storing  the  information 
is  the  same  as  described  for  LBLX.  The  maximum 
label  length  is  24  characters. 
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Restrictions: 

1. 


The  value  of  XHIGH  must  be  equal  to  or  greater 
than  XLOW. 

The  value  of  YHICH  must  bo.  equal  to  or  greater  than  YLOW. 

Labelling  information  will  be  truncated  if  the  maximum 
label  length  is  exceeded;  i.e.,  30  characters  cr  x an 

24  for  y. 
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PPLOT 


Subroutine  Specification: 

SUBROUTINE  PPLOT  (IPLOT,  ICHR,  NPTS , XV,  YV) 
LOGICAL  IPLOT 
INTEGER  ICHR,  NPTS 
DIMENSION  XV(1),  YV(1) 


Purpose: 

To  store  plot  data  and  initiate  plotting.  The  IPLOT 

parameter  controls  plot  initiation. 

Usage: 

The  Fortran  calling  sequence  is: 

CALL  PPLOT  (IPLOT,  ICHR,  NPTS,  XV,  YV'' 

where : 

IPLOT  ■ LOGICAL  control  variable  for  plotting.  Specify 

•FALSE,  if  d«ta  values  are  to  be  stored  and  plot- 
ting is  to  be  deferred;  (used  for  multiple  plots  per 
graph).  Specify  .TRUE,  if  data  values  are  to  be  stored 
and  plotted. 

ICHR  ■ a single  Hollerith  character  to  be  used  as  the  plot 
symbol;  (e.g.,  1H*) . 

NPTS  • number  of  data  points  to  be  plotted. 

XV, YV  ■ vectors  containing  coordinate  values  (x  ,y.)  to  be 
plotted. 

Restrictions : 

1.  If  more  than  one  data  point  occupies  the  same  plot 
position,  the  most  recent  value  processed  will  be  used. 

2.  Data  values  outside  the  range  specified  by  PSETUP  will 
not  be  plotted.  They  will  be  listed  below  the  plot. 
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SECOND 

Subroutine  Specification: 

SUBROUTINE  SECOND(ISEC) 

INTEGER  I SEC 

Programmed  by:  L.  C.  Lintner 

L!S2S&''  This  routine  provides  the  user  with  a means  of  determining 
the  amount  of  central  Processor  (CPU)  time  that  he  has  used  Th  time 
is  accumulated  from  the  beginning  of  the  run;  i.e..  it  includes  the 
time  accrued  for  all  preceding  tasks  in  the  run  setup. 


Usage: 


SECOND  is  an  assembly  language  routine  in  "FORTRAN  callable"  form. 


The  FORTRAN  calling  sequence  is: 

CALL  SECOND (I SEC) 

The  parameter  ISEC  must  be  an  integer  variable.  The  result 
returned  in  ISEC  is  the  integer  number  of  CPU  seconds  used  since  the  star. 

of  the  run. 

This  routine  may  also  be  called  as  an  integer  function: 

L » NSEC (ISEC) 

In  this  case,  the  integer  number  of  seconds  will  be  returned  as  the  value 
of  the  function  in  addition  to  being  stored  in  ISEC. 
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